“Orden de red, morfometría y análisis hortoniano usando r.stream*"

#

Imprimir lista de mapas ráster y vectoriales dentro en la región/localización activa

• Nótese que los paquetes requeridos en esta sessión (rgrass7, raster, leaflet, leafem), fueron en el bloque anterior al ejecutarse el código contenido en el archivo extraer-red-de-drenaje-con-r-stream.Rmd. Igualmente, dicho bloque de código creó todos los objetos necesarios para realizar este tutorial.

execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/LfpNetwork-flds-Oz
## raster/LfpNetwork-outlet-Oz
## raster/MASK
## raster/accum-de-rwshed
## raster/aspect
## raster/basins
## raster/dem
## raster/drainage-dir-de-rstr
## raster/drainage-dir-de-rwshed
## raster/half-basins
## raster/order-hack-gravelius
## raster/order-horton
## raster/order-shreve
## raster/order-strahler
## raster/order-topology
## raster/ozama-basin
## raster/ozama-stream-de-rstr
## raster/pcurv
## raster/r-stream-basins-1
## raster/r-stream-basins-2
## raster/r-stream-basins-3
## raster/r-stream-basins-4
## raster/r-stream-basins-5
## raster/r-stream-basins-6
## raster/r-stream-basins-7
## raster/slope
## raster/stream-de-rwshed
## raster/tcurv
## vector/LfpNetwork_lfp_all_final_ozm
## vector/LfpNetwork_lfp_ozm
## vector/LfpNetwork_outlet_Oz
## vector/LfpNetwork_tributaries_ozm
## vector/LfpNetwork_tributaries_preconf_ozm
## vector/c_ozama
## vector/dem_extent
## vector/order_all
## vector/ozama_basin
## vector/ozama_stream_de_rstr
## vector/r_stream_basins_1
## vector/r_stream_basins_2
## vector/r_stream_basins_3
## vector/r_stream_basins_4
## vector/r_stream_basins_5
## vector/r_stream_basins_6
## vector/r_stream_basins_7

Crear mapa de dirección de flujo a partir de r.stream

execGRASS(
  "r.stream.extract",
  flags = c('overwrite','quiet'),
  parameters = list(
    elevation = 'dem',
    threshold = 80,
    direction = 'drainage-dir-de-rstr'
  )
)

Crear mapas de órdenes de red

execGRASS(
  "r.stream.order",
  flags = c('overwrite','quiet'),
  parameters = list(
    stream_rast = 'ozama-stream-de-rstr',
    direction = 'drainage-dir-de-rstr',
    elevation = 'dem',
    accumulation = 'accum-de-rwshed',
    stream_vect = 'order_all',
    strahler = 'order-strahler',
    horton = 'order-horton',
    shreve = 'order-shreve',
    hack = 'order-hack-gravelius',
    topo = 'order-topology'
  )
)
## Warning in execGRASS("r.stream.order", flags = c("overwrite", "quiet"), : The command:
## r.stream.order --overwrite --quiet stream_rast=ozama-stream-de-rstr direction=drainage-dir-de-rstr elevation=dem accumulation=accum-de-rwshed stream_vect=order_all strahler=order-strahler horton=order-horton shreve=order-shreve hack=order-hack-gravelius topo=order-topology
## produced at least one warning during execution:
## WARNING: Vector map <order_all> already exists and will be overwritten
## WARNING: Vector map <order_all> already exists and will be overwritten

Mostrar lista nuevamente

execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/LfpNetwork-flds-Oz
## raster/LfpNetwork-outlet-Oz
## raster/MASK
## raster/accum-de-rwshed
## raster/aspect
## raster/basins
## raster/dem
## raster/drainage-dir-de-rstr
## raster/drainage-dir-de-rwshed
## raster/half-basins
## raster/order-hack-gravelius
## raster/order-horton
## raster/order-shreve
## raster/order-strahler
## raster/order-topology
## raster/ozama-basin
## raster/ozama-stream-de-rstr
## raster/pcurv
## raster/r-stream-basins-1
## raster/r-stream-basins-2
## raster/r-stream-basins-3
## raster/r-stream-basins-4
## raster/r-stream-basins-5
## raster/r-stream-basins-6
## raster/r-stream-basins-7
## raster/slope
## raster/stream-de-rwshed
## raster/tcurv
## vector/LfpNetwork_lfp_all_final_ozm
## vector/LfpNetwork_lfp_ozm
## vector/LfpNetwork_outlet_Oz
## vector/LfpNetwork_tributaries_ozm
## vector/LfpNetwork_tributaries_preconf_ozm
## vector/c_ozama
## vector/dem_extent
## vector/order_all
## vector/ozama_basin
## vector/ozama_stream_de_rstr
## vector/r_stream_basins_1
## vector/r_stream_basins_2
## vector/r_stream_basins_3
## vector/r_stream_basins_4
## vector/r_stream_basins_5
## vector/r_stream_basins_6
## vector/r_stream_basins_7

Visualizar la red con leaflet

Simbología única

order <- readVECT('order_all')

order4326 <- spTransform(order, CRSobj = CRS("+init=epsg:4326"))
leaflet() %>% 
  addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
  addPolylines(
    data = order4326, weight = 3, opacity = 0.7, group = 'order',
    label = ~as.character(strahler),
    highlightOptions = highlightOptions(color = "white",
                                      weight = 5, bringToFront = F, opacity = 1),
    labelOptions = labelOptions(noHide = F,
                                style = list(
                                  "font-size" = "8px",
                                  "background" = "rgba(255, 255, 255, 0.5)",
                                  "background-clip" = "padding-box",
                                  "padding" = "1px"))) %>% 
  leafem::addHomeButton(extent(order4326), 'Ver todo') %>% 
  addLayersControl(
    overlayGroups = c('terrain','order'),
    options = layersControlOptions(collapsed=FALSE))

Simbología aplicando grosor según orden de red

orden_de_red <- leaflet() %>% 
  addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
  addPolylines(
    data = order4326, weight = order4326$strahler*1.5, opacity = 0.7, group = 'order',
    label = ~as.character(strahler),
    highlightOptions = highlightOptions(color = "white",
                                      weight = 5, bringToFront = F, opacity = 1),
    labelOptions = labelOptions(noHide = F)) %>% 
  leafem::addHomeButton(extent(order4326), 'Ver todo') %>% 
  addLayersControl(
    overlayGroups = c('terrain','order'),
    options = layersControlOptions(collapsed=FALSE))
orden_de_red

orden_de_red %>% mapview::mapshot(file = 'orden_de_red_salida.png')

Delimitar cuencas según orden de red de Strahler

Obtener órdenes de red mínimo y máximo

#Estadísticas para obtener los valores mínimo y máximo del orden de red de Strahler
rinfo.ordstra <- execGRASS(
  'r.info',
  flags = 'r',
  parameters = list(
    map = 'order-strahler'
  )
)
## min=1
## max=7
#Órdenes de red mínimo y máximo
minmaxord <- as.numeric(
  stringr::str_extract_all(
    attributes(rinfo.ordstra)$resOut,
    "[0-9]+"
  )
)
minmaxord
## [1] 1 7

Delimitar cuencas, convertirlas de ráster a vectorial

sapply(
  min(minmaxord):max(minmaxord),
  function(x){
    execGRASS(
      "r.stream.basins",
      flags = c('overwrite','c','quiet'),
      parameters = list(
        direction = 'drainage-dir-de-rstr',
        stream_rast = 'order-strahler',
        cats = as.character(x),
        basins = paste0('r-stream-basins-',x)
      )
    )
    execGRASS(
      "r.to.vect",
      flags=c('overwrite','quiet'),
      parameters = list(
         input = paste0('r-stream-basins-',x),
         output = paste0('r_stream_basins_',x),
         type = 'area'
      )
    )
  }
)
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-1 output=r_stream_basins_1 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_1> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_1> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-2 output=r_stream_basins_2 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_2> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_2> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-3 output=r_stream_basins_3 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_3> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_3> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-4 output=r_stream_basins_4 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_4> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_4> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-5 output=r_stream_basins_5 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_5> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_5> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-6 output=r_stream_basins_6 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_6> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_6> already exists and will be
##          overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-7 output=r_stream_basins_7 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_7> already exists and will be
##          overwritten
## WARNING: Vector map <r_stream_basins_7> already exists and will be
##          overwritten
## [1] 0 0 0 0 0 0 0

Representar las cuencas con leaflet

sapply(
  min(minmaxord):max(minmaxord),
  function(x){
    assign(
      paste0('orden', x),
      spTransform(readVECT(paste0('r_stream_basins_',x),driver = 'SQLite'), CRSobj = CRS("+init=epsg:4326")),
      envir = .GlobalEnv)
  }
)

paleta <- RColorBrewer::brewer.pal(12, 'Set3')
cuencas_y_orden_de_red <- leaflet() %>% 
  addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
  addPolygons(data = orden7, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O7') %>%
  addPolygons(data = orden6, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O6') %>%
  addPolygons(data = orden5, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O5') %>%
  addPolygons(data = orden4, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O4') %>% 
  addPolygons(data = orden3, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O3') %>%
  addPolygons(data = orden2, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O2') %>%
  addPolygons(data = orden1, stroke = T, weight = 2,
              color = ~paleta, fillOpacity = 0.4, group = 'O1') %>%
  addPolylines(
    data = order4326, weight = order4326$strahler*1.5,
    opacity = 0.7, group = 'str_order') %>%
  leafem::addHomeButton(extent(order4326), 'Ver todo') %>% 
  addLayersControl(
    overlayGroups = c('terrain','O1','O2','O3','O4','05','06','07','str_order'),
    options = layersControlOptions(collapsed=FALSE))
cuencas_y_orden_de_red

cuencas_y_orden_de_red %>% mapview::mapshot(file = 'cuencas_y_orden_de_red_salida.png')

Estadísticas de red resumidas por orden de red.

execGRASS(
  "r.stream.stats",
  flags = c('overwrite','quiet','o'),
  parameters = list(
    stream_rast = 'order-strahler',
    direction = 'drainage-dir-de-rstr',
    elevation = 'dem',
    output = 'ozama_stats.txt'
  )
)
file.show('ozama_stats.txt')
d <- read.csv("ozama_stats.txt", skip=1, header=TRUE)
d
##   order num_of_streams avg_length    avg_area avg_slope avg_grad
## 1     1           1328   1.157222    1.531833  0.019501 0.013230
## 2     2            274   2.787560    6.990865  0.013704 0.007789
## 3     3             63   8.174857   36.891303  0.009048 0.004075
## 4     4             13  18.228154  180.416363  0.007264 0.001903
## 5     5              5  21.793405  494.576211  0.004899 0.000828
## 6     6              2  17.109400 1364.385565  0.002296 0.000297
## 7     7              1  23.929262 3376.766977  0.005338 0.000190
##   avg_elev.diff sum_length sum_area
## 1     15.912760 1536.79041 2034.274
## 2     21.873087  763.79156 1915.497
## 3     39.441891  515.01597 2324.152
## 4     45.280315  236.96600 2345.413
## 5     12.680176  108.96702 2472.881
## 6      5.009900   34.21880 2728.771
## 7      4.555471   23.92926 3376.767
plot(num_of_streams~order, data=d, log="y")
mod <- lm(log10(num_of_streams)~order, data=d)
abline(mod)
text(2, 20, 'logN=2.064-0.544u')

rb <- 1/10^mod$coefficients[[2]]
rb
## [1] 3.361632

Estadísticas de red ampliadas

execGRASS(
  "r.stream.stats",
  flags = c('overwrite','quiet'),
  parameters = list(
    stream_rast = 'order-strahler',
    direction = 'drainage-dir-de-rstr',
    elevation = 'dem',
    output = 'ozama_stats_expanded.txt'
  )
)
file.show('ozama_stats_expanded.txt')

Limpiar archivo de bloqueo del conjunto de mapas de GRASS

rbp <- mean(d$num_of_streams[-length(d$num_of_streams)]/d$num_of_streams[-1])
rbp
## [1] 3.523679

unlink_.gislock()